Mon. Not. R. Astron. Soc. OOO.fTHTSI (2011) Printed 27 January 2013 (MN JATfcX style file v2.2) 



Rapid inward migration of planets formed by 
gravitational instability 



Clement Baruteau 1 ' 2 * , Farzana Mem 3,4 and Sijme-Jan Paardekooper 1 

1 DAMTP, University of Cambridge, Wilberforce Road, Cambridge CBS OWA, UK 

2 Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA 

^School of Physics, University of Exeter, Stacker Road, Exeter, EX4 4QL, UK 

4 Institut fur Astronomie und Astrophysik, Universitat Tubingen, Auf der Morgenstelle 10, 72076 Tubingen, Germany 



o 

(N 
G 

(N 

w 

6 



> 

00 

^- 
o 

\6 
o 



Accepted 2011 June 1. Received 2011 May 30; in original form 2011 March 25 



ABSTRACT 

The observation of massive exoplanets at large separation (> 10 AU) from 
their host star, like in the HR 8799 system, challenges theories of planet for- 
mation. A possible formation mechanism involves the fragmentation of massive 
self-gravitating discs into clumps. While the conditions for fragmentation have 
been extensively studied, little is known of the subsequent evolution of these 
giant planet embryos, in particular their expected orbital migration. Assuming 
a single planet has formed by fragmentation, we investigate its interaction with 
the gravitoturbulent disc it is embedded in. Two-dimensional hydrodynamical 
simulations are used with a simple prescription for the disc cooling. A steady 
gravitoturbulent disc is first set up, after which simulations are restarted includ- 
ing a planet with a range of masses approximately equal to the clump's initial 
mass expected in fragmenting discs. Planets rapidly migrate inwards, despite the 
stochastic kicks due to the turbulent density fluctuations. We show that the mi- 
gration timescale is essentially that of type I migration, with the planets having 
no time to open a gap. In discs with aspect ratio ~ 0.1 at their forming loca- 
tion, planets with a mass comparable to, or larger than Jupiter's can migrate 
in as short as 10 4 years, that is, about 10 orbits at 100 AU. Massive planets 
formed at large separation from their star by gravitational instability are thus 
unlikely to stay in place, and should rapidly migrate towards the inner parts of 
protoplanetary discs, regardless of the planet mass. 

Key words: accretion, accretion discs — turbulence — methods: numerical — 
planetary systems: formation — planetary systems: protoplanetary discs 
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1 INTRODUCTION 

The large majority of exoplanets have been discovered 
by the radial velocity and transiting techniques. Re- 
cent progress in the adaptive optics technique has re- 
vealed about 20 exoplanets by direct imaging. Since the 
latter requires a high brightness contrast between the 
planet and the parent star, massive planets have been 



observed at separations >. 10 AU (e.g.. | Kalas et al.ll2008l; 
Marois et al.ll2008l;lLagrange et aI]|201oT Lafreniere et al.l 



20101 : 1 Marois et a l. 2010). One s uch example is the HR 
8799 system (|Marois et alj|20ich . which comprises four 
planets with masses evaluated between 7 and 10 Jupiter- 
mass, and estimated separations of 14, 24, 38 and 68 AU, 
respectively. 
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How can the large orbital separations of such mas- 
sive planets be explained? Wit h the standard core ac- 
cretion formation scenario (e.g., Safronov| | l969l : iMizunol 
Il98d : iPollack et ail 1 19961 : lldafc Linl |2004| ). it is diffi- 
cult to form Jupiter-like planets in isolation further 
than ~ 10 AU from a Sun-like star. Still, the forma- 
tion of a first ga p-opening planet close to the snow line 
(|lda fc Linll2004l ) could pave the way for the formation 
of a second-genera tion planet at the o uter edge of the 
first planet's ga p (jBrvden et all l200(j iThommed 120051 : 
llda fc Lin| [2008). Provided that the second planet opens 
a gap, its outer edge could also trigger the formation of 
a third-generation planet. Subsequent planet formation 
could proceed as a sequence, as long as there are enough 
planetesimals and gas left. 

The tidal interaction between a planet and 
its nascent protoplanetary disc alters the planet's 



© 2011 RAS 



2 C. Baruteau et al. 



semi-major axis, causing its orbital migration 
dGoldreich fc Tremaind [1980). Assuming the central 
object is a Sun-like star, a planet typically more massive 
than Saturn opens a gap around its orbit and generally 
migrates inwards. It is thus unlikely that a massive 
planet formed through the core-accretion scenario within 
~ 10 AU of its host star could migrate to a separation as 
large as 100 AU. A notable exception was proposed by 
ICrida et al. (2009), who showed that the tidal interaction 
between a pair of resonant massive planets embedded 
in a common gap, and their protoplanetary disc, could 
drive both planets to orbital separations much larger 
than that of their birth place, providing some specific 
conditions on the planets mass ratio and on the disc's 
turbulent viscosity. 

Another in-situ formati on scenario is based on th e 
gravitational instability (e.g.. lCamero ni l 19781 : iBossll 19971 ). 
It involves early-stage protoplanetary discs, which are 
massive enough for its self-gravity to play a part in its 
evolution. The probability that discs form fragments is 
described by two dimensionless parameters. One is the 
Toomre Q-parameter, Q = c s K/7rGE, where G is the 
gravitational constant, c a and E denote the disc's sound 
speed and surface mass density, and k is the epicyclic fre- 
quency. For Kepleria n discs, k appr oximately equals the 
angular frequency f2. iToonirel (1964) showed that, for an 
innnitesimally thin disc to fragment, Q < 1. The second 
parameter describing the stability of a self-gravitating 
disc is the disc's cooling timescale in units of the orbital 
timescale, /3 = t coo \ Q, where t coo i = e(de/dt)~ 1 is the 
cooling timesca le and e the disc's thermal energy density 
l|Gammiell200ll ). The critical value of the cooling param- 
eter /3, below which fragmentation wil l occur if Q < 1 , 
depends on the disc's adiabatic index (jRice et al] 20051). 
and on the local disc to primary mass ratio ( Meru & Bate 
l2011bl). Its exac t value is currently somewhat uncertain 
as lMeru fc Bate! (|201 lal ) have recently shown that previ- 
ous simulations that determined this critical value were 
unresolved. However, current simulations show that the 
critical value for could be in the range 3 — 18 (jGammiel 
12001] ; iRice et al]|2005l : iMeru fc Batdl2011al lbh. 

The vast majority of studies of planet formation by 
gravitational instability have focused on the conditions 
to fragment discs into bound objects. It is thought to op- 
erate at separations larg er than 30 — 50 AU from the cen- 
tral (Sun-like) star (e.g..lRafikoyll2005l : lMat zner fc Levinl 
120051 ; IStamatellos fc Whitworthll2008l ). However, little is 
known of the evolution of these objects, in particular 
their expected migration resulting from the gravitational 
interaction with the gravitoturbulent disc they are em- 
bedded in. A few numerical stu dies have observed that 
clumps could drift inwar ds ie.g., iMaver et al1l2002l : iBossI 
120051 ; ICha fc NavakshinllioibT) . and some of them have 
focused on the early phases of disc formation and evolu- 
tion following the collap s e from th e prestellar core stag e 
l|Vorobvov fc Basul l2006t l2010al lbt iMachida et ail l201ll ). 
As these giant planet embryos migrate inwards, they 
may become tidally disrupted (when their radius becomes 
comparable to the clump's Hill radius), thereby delivering 
a vari ety of planets in the inner parts of protoplanetary 
discs l)Navakshinll2010h . 

In this paper, we examine the orbital evolution of 



a single planet embedded in its nascent gravitoturbulent 
disc, following the assumption that the planet has formed 
by gravitational instability. We show that the planet mi- 
grates inw ards in a timescale very similar to th at of type I 
migration ( Paardcko oper fc Papaloizoul f2009) in the ab- 
sence of turbulence. Our physical and numerical models 
are described in Section [2] In Section we present our 
results of simulations. Conclusions and future directions 
are drawn in Section [4] 



2 MODEL 

We investigate the tidal interaction between a massive 
planet and the protoplanetary disc it formed in through 
the gravitational instability scenario. Once a gaseous 
clump has formed by fragmentation, it contracts until 
H2 dissociation occurs, causing a rapid collapse of the 
clump down to planetary sizes. Our study does not ad- 
dress the formation of the planet, and we consider the 
clump and the planet as one object. In the following, our 
system of interest (clump and planet) is referred to as the 
planet. Our study focuses on its migration driven by the 
interaction with the disc, following the assumption that 
the planet has formed by gravitational instability. As de- 
scribed below, our approach is to first set up a steady 
state gravitoturbulent disc, then include a single planet 
and follow its orbital migration. 

The disc's properties after fragmentation are poorly 
constrained. They depend on how many clumps form, 
and how massive they are. In our study, we assume that 
the disc properties (density, temperature) are virtually 
unchanged by the fragmentation stage. In particular, we 
assume that the disc remains gravitoturbulent. This as- 
sumption is motivated by the findings of previous nu- 
merica l studies on disc fragmentation (e.g.. lMeru fc Bate] 
l2011bl ). where the images of the discs show similar struc- 
tures after the formation of a clump. 

We carried out two-dimensional (2D) hydrodynam- 
ical simulations using the g rid-based, staggered-mesh 
code FARGC0 (M assetl [2000). for which a self-gravity 
solver based on fast Fourier transforms was presented in 
iBaruteau fc Massed (|2008bl) . When calculating the self- 
gravitating acceleration, all mass is assumed to be con- 
fined to the disc midplane (the razor-thin disc approxima- 
tion) . The thermal energy density e satisfies the equation 



(ev) 



-(7 - l)eV • v + 1 



"7~cool 



(i) 



where v denotes the gas velocity, and 7 is the adi- 
abatic index (its value can b e mapped to a three- 
dimensional adiabatic index, see lGammidl200ll . and ref- 
erences therein) . In Eq. ([1} , the cooling timescale r coo i = 
/3f2 _1 , with /3 taken to be constant, and the heating 
source term Q^ ulk is provided by artificial viscous heat- 
ing at shocks arising from the gravitoturbulence. A von 
Neumann-Richtmyer artificial bulk viscosity is used, as 
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Figure 1. Azimuthally-averaged surface density (left panel), temperature (middle panel), and Toomre Q-parameter (right panel) 
obtained at 30 orbits for our three values of = t coo \ Q (solid curves). The azimuthal average of Q was calculated as (c s ) (Q)/irG{T,), 
where (.) stands for the azimuthal average. Dashed curves depict the initial value of each quantity. 



described in IStone fc" Normanl {l992), where the coeffi- 
cient C2 is taken equal to 1.4 (C2 measures the num- 
ber of zones over which a shock is spread over by the 
artificial viscosity). Providing angular momentum trans- 
port driven by gravitoturb ulence occurs locally, which 
is not necessarily the case <|Balbus fc Papaloizoul [l999; 
ICossins et al.|[2009l ). the alpha parameter associ ated with 
the a ngular momentum flux density is given by l|Gammiel 

mi) 



(2) 



9 7(7-1)/?' 



which indicates that varying ft is equivalent to varying the 
viscous stress in the disc. Preliminary calculations with 
our fiducial setup (initial profiles and grid's resolution are 
described below) showed that the critical value of /3 below 
which fragmentation occurs ranges from 10 to 15. Since 
we require the disc to be in a gravitoturbulent state, but 
we do not require further fragmentation so as to focus on 
the migration of a single planet, we simulated three disc 
models with /3 = 15, 20 and 30. 

Each of our three models initially comprises a O.4A/0 
2D disc around a star of fixed mass = IMq, and 
spanning a radial range 20 < r < 250 AU. (As will be 
shown in Section [3.11 the disc mass rapidly decreases to 
about 0.2 - 0.25A'/ Q after 30 orbits at 100 AU). The ini- 
tial surface mass density decreases as r~ 2 , and it equals 
~ 20 g cm -2 at 100 AU. The initial disc's temperature 
is ~ 15 K at 100 AU, and it decreases as r^ 1 (we as- 
sume a mean-molecular weight fi = 2.4 and an adiabatic 
index 7 = 5/3). This corresponds to setting the initial 
disc aspect ratio h (pressure scale height to radius ratio 
H/r) equal to 0.1 at 100 AU. The disc's initial Toomre 
Q-parameter is thus uniform, equal to 2. A small level 
(~ 0.1%) of white noise is added initially to break the 
disc's axisymmetry. Here and in the following, disc and 
planet quantities are expressed in physical units. Since 
our simulations are scale- free, all the results presented 
throughout this paper can be easily rescaled by using 
different units of mass, length and temperature. 

In all our simulations, the grid is covered with 512 
and 1536 cells along the radial and azimuthal directions, 
respectively. A logarithmic spacing is used along the ra- 
dial direction, as required by the self-gravity module in 



FARGO. Grid cells are approximately square, and the 
initial pressure scale height is resolved by about 20 cells 
along each direction. The radial and azimuthal compo- 
nents of the self-gravitating acceleration are smoothed 
over a softening length, e sg . We take e sg = 3 x 10~ 4 r 
throughout this study, which is a factor 3 x 10 -3 times 
the initial pressure scale height. The influence of a large 



of finite thickness, will be examined in Sect. I3.2~4l where 
we find that the choice for e sg does not significantly alter 
the orbital evolution of a planet embedded in a gravito- 
turbulent disc. Outflow boundary conditions are used at 
the grid's inner and outer edges. Hydrodynamical equa- 
tions are solved in the frame centred onto the central star. 
Since a non-inertial frame of reference is adopted, the in- 
direct terms in the expression for the disc's gravitational 
potential are included. 



3 RESULTS OF HYDRODYNAMICAL 
SIMULATIONS 

Our numerical simulations were carried out in two steps. 
In a first step, a gravitoturbulent disc model was set 
up, whose properties are described in Section 13.11 We 
then restarted the simulations, including planets of sev- 
eral masses, and following the time-evolution of their or- 
bital separation. The results of the restart simulations 
are presented in Section 13.21 Further investigation on the 
impact of the disc's temperature follows in Section [3.31 



3.1 Gravitoturbulent stage 

We describe in this section the results of our three 
disc models with cooling parameters /3 — 15, 20 and 
30. The absence of heating source initially leads to de- 
crease the disc's temperature, and thus the Toomre Q- 
parameter, until gravitoturbulence sets in. As both the 
initial Q-profile and /3 are taken to be uniform (the 
cooling timescale thus scales proportional to the or- 
bital timescale), gravitoturbulence progressively devel- 
ops through the disc, starting from the inner edge. After 
a few orbits at the disc's outer edge, a thermal quasi- 
equilibrium is reached over the whole disc, where cooling 
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Figure 2. Gravitational, Reynolds and total alpha parame- 
ters, azimuthally- and time-averaged from 30 to 60 orbits at 
100 AU (that is, over S3 8 orbits at 250 AU, the location of 
the disc's outer edge). Results are displayed for the disc model 
with cooling parameter /3 = 20. The dashed curve shows the 
total alpha parameter expected from viscous transport theory, 
given by Eq. f2l . 



is approximately balanced by shock heating due to grav- 
itoturbulence. 

The disc's initial surface density that we assume 
does not correspond to a mechanical equilibrium with 
the gravitoturbulent stress. The disc therefore adjusts 
both its surface density and temperature profiles to 
satisfy a mechanical equilibrium along with a uniform 
Toomre-Q profile. This is illustrated in Fig. [T] where the 
azimuthally-averaged surface mass density, temperature, 
and Toomre Q-parameter are depicted at 30 orbits at 100 
AU. Their initial value is shown with a dashed curve. 
In this quasi-steady state, the surface density profile de- 
creases approximately as r~ 3 ^ 2 , and the disc's temper- 
ature is almost uniform, as is the Toomre Q-parameter. 
We comment that the high initial disc mass (OAMq) pro- 
motes global modes that, along with the outflow bound- 
ary condition at the grid's inner edge, lead to a rapid 
decrease in the disc mass. After 30 orbits at 100 AU, the 
latter ranges from 0.21 Mo to 0.24M© for /3 increasing 
from 15 to 30. The rate at which the disc mass decreases 
is maximum at the onset of gravitoturbulence, and then 
declines (but does not vanish) as a quasi-steady state is 
attained. 

The alpha parameters associated with the Reynolds 
and gravitational turbulent stresses are denoted by ctRey 
and a sg , respectively. Details of their calculation in a 2D 
disc model are given in the Appendix. Fig. [2]displays the 
azimuthally- and time-averaged radial profiles of QRcy, 
a sg and of atot = ct ss + ciRcy • Time average is done from 
30 to 60 orbits at 100 AU (that is, over « 8 orbits at the 
grid's outer edge). Results are shown for j3 = 20. We note 
that Q sg and «R e y take very similar value s over the whole 
disc. The same behavior was obtained bv lGammiel l|200ll) 
with 2D local shearing-sheet calculations. The decrease 
in both a sg and an ay at large (> 200 AU) separation 
is most probably due to the decreasing grid's resolution, 
as a logarithmic spacing is used along the radial direc- 
tion. We comment that a to t is locally up to 50% larger 



than the value expected from viscous (local) transport 
theory, given by Eq. The difference is probably due 
to our large disc mass (> 0.2M*) and large aspect ratio 
(h ~ 0.1), which ten d to foster global tran sport of angu- 
lar momentum (e.g.. lLodato fc Ricell2004h . Although not 
shown here, very similar results were obtained for ft = 15 
and P = 30. 

Finite resolution effects (numerical diffusion) leads 
to a negligible transport of angu lar momentum compared 
to gravitoturbulence (see e.g., iLin fc Papaloizoul |2010| , 
who estimated the numerical "alpha viscosity" in the 
FARGO code to be < 10~ 5 for a number of grid cells sim- 
ilar to ours). Artificial bulk viscosity most likely leads to 
a small, if not negligible transport of angular momentum 
compared to gravitoturbulence, but note that it is unclear 
how its contribution can be isolated in our simulations. 
We point out that if artificial bulk viscosity induced a 
net transport of angular momentum on average, it would 
lead to a discrepancy between the value of a sg + QR cy 
measured in local, steady-s tate disc models and the al- 
pha value given in Eq. p. lGammiel l|200ll ) carried out 
2D local shearing-sheet simulations of gravitoturbulent 
disc models using the ZEUS code, which is very simi- 
lar to the FARGO code. In particular, both codes use 
the same prescription for the artificial bulk viscosity. He 
found that a sg + evR oy agrees with the value given by the 
expression in Eq. ([2} to within ~ 10%, which shows that 
artificial bulk viscosity should lead to a small transport of 
angular momentum compared to gravitoturbulence. Also, 
although not shown here, we have checked that increas- 
ing t he amount of artificial b ulk viscosity (the coefficient 
C2 in lStone fc Normanlll992l ) by a factor of 3 leads to a 
negligible change in the values of ttR cy and a sg . 



3.2 Restart with an embedded planet 



We restarted the simulations of Section [3.11 at 30 orbits, 
including a single planet at 100 AU. The disc's tem- 
perature and surface mass density at this location are 
close to their initial value for our three values of /3 (see 
Fig. [TJ. Note in particular that the disc aspect ratio at 
the planet's initial location is ~ 0.1. Akin to the clump's 
contraction timescale, the final mass of pla nets formed 
by gravitational instability is uncert ain (e.g.. lBolev et all 
l2010l ; iHelled fc Bodenheimerl l201ll ), and may substan- 
tially differ f rom the clump's in itial mass. The latter is 
estimated bv lBolev et all l|2010l ) as ~ h% M* for Q ~ 1.5, 
with hp the disc aspect ratio at the clump's forming lo- 
cation. This estimate of the clump's initial mass, which 
is based on the fragmentation of spiral arms, is typi- 
cally a factor 5 — 10 smaller than the Toomre mass (see 
iBolev et al]|2010l . their figure 1, and references therein). 
We considered three planet masses corresponding to a 



2 Similar gravitoturbulent disc models, but with reduced 
aspect ratio in a s teady state, have been simulated by 



iPaardekooper et al. (2011) with the same hydrodynamical 
code and similar grid resolution. Using the same method to 
calculate the stresses as described in the appendix, they find 
an averaged value of atot in very good agreement with the 
value given in Eq. J2j. 
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Figure 3. Contours of the disc's surface density obtained in 
a restart simulation with = 30, where a Jupiter-mass planet 
has been introduced at 100 AU. Results are displayed three 
orbits after the restart time. The planet is now located at 
about 60 AU from the central star (at x ~ —60 AU, y ~ 60 
AU). 

planet to primary mass ratio q = M p /M+ around /i 3 : a 
Saturn-mass planet (q — 3 x 10 -4 , q/hp ~ 0.3), a Jupiter- 
mass planet (q — 10~ 3 , q/hp ~ 1), and a 5 Jupiter-mass 
planet (q = 5 x 10 -3 , q/hp ~ 5). The planets gravi- 
tational potential is softened over a smoothing length 
e — 3 x 10~ 2 r p , with r p the planet's orbital separation, 
and it takes its maximum value at the restart time (we 
do not gradually ramp up the planet's mass). We will 
show in Section 13.2.41 that our results are not altered 
by a more gentle introduction of the planet's potential 
in the disc over a dynamical timescale, the typical for- 
mation timescale of a clump. We comment that in self- 
gravitating discs, the mass that is relevant for disc-planet 
interactions is the mass of the planet and of the gas enve- 
lope surrounding it. Although for simplicity the planet's 
mass is fixed, gas can be accreted onto its envelope (cir- 
cumplanetary disc), thereby contributing to increasing 
the effective planet's mass. Note however that this ac- 
cretion rate is probably not realistic due to our simple 
prescription for the disc cooling. Since the disc is fully 
self-gravitating, the calculation of the force exerted on the 
planet takes all the disc into ac count, it does not partly 
exclude the planet's Hill radius ()Crida et al.ll200^ ). 

Fig. [3] displays the surface density obtained for the 
disc model with ft = 30, perturbed by a Jupiter-mass 
planet. Density contours are shown only three orbits af- 
ter the restart, the planet is now located at « 60 AU from 
the central object. The density perturbation due to the 
planet's potential is comparable to the turbulent density 
perturbations. This similarity stresses that the magni- 
tude of the stochastic torque acting on the planet can be 
similar to that of the disc-planet tidal torque obtained 
in the absence of turbulence (e.g., [Nelson fe Papaloizoul 



|2004 iBaruteau fc"Lir]|2010h . In such turbulent configu- 
rations, the migration timescale must be evaluated in a 
statistical way, by running several simulations. For each 
planet mass and each value of /3, we performed a suite of 
eight simulations, varying the planet's azimuth at restart 
as Lp v = 7r/4 x i, i £ [0 — 7] (that is, independently of 
the disc's density structure at restart). The migration 
timescales obtained for these series of runs are presented 
in Section ^. 2. II The impact of stochastic density fluctua- 
tions is described in Section ^. 2. 21 A comparison with pre- 
dictions of laminar disc models follows in Section 13.2.31 
and a brief discussion on the impact of numerics is given 
in Section WTM 



3.2.1 Migration timescale 

The time evolution of the planets' orbital separation is 
displayed in Fig. [4] Planet mass increases from left to 
right, and the cooling parameter f3 increases from top to 
bottom. In all panels, time is expressed in years (bottom 
x-axis), and in orbital periods at 100 AU (top x-axis). Re- 
sults are shown for the eight evenly-spaced values of the 
planet's azimuth ip p at restart (increases with increasing 
color's wavelength). 

The net trend coming out of our results is that, in 
spite of the stochastic kicks triggered by turbulence, plan- 
ets migrate inwards very rapidly. Averaging over (p v and 
P, the migration timescale (defined as the time to migrate 
from 100 AU to 20 AU, the location of the grid's inner 
edge) is typically as short as 3 orbits for the 5 Jupiter- 
mass planet, 6 orbits for the Jupiter-mass planet, and 
about 8 orbits for the Saturn-mass planet. In addition, at 
fixed P, the amplitude of stochastic kicks increases with 
decreasing planet mass, and so does the spread in the 
migration timescale. Both results may be qualitatively 
interpreted as follows. The torque exerted by the disc 
on the planet may be decomposed as a background tidal 
torque, plus a stochastic torque due to turbulent fluctu- 
ations. In the mass regime that we consider (q/hp ~ 1), 
the amplitude of the tidal torque should increase with 
increasing planet mass. However, like in discs with tur- 
bulence driven by the magneto-rotational instability, it is 
unclear if the value of the tidal torque is similar to that 
predi cted in laminar disc models, in the absence of turbu- 
lence (|Nelson fc Papa loizou 200J) . A detailed comparison 
with results of laminar disc models will be presented in 
Section 13.2.31 The standard deviation of the stochastic 
torque's distribution is primarily controlled by atot, and 
therefore by the cooling parameter /3 (see Eq. [2j. For a 
given /3, the smaller the planet mass, the smaller the tidal 
torque to stochastic torque ratio, the more sensitive the 
planet is to stochastic kicks, and therefore the larger the 
spread in the migration timescale. In the same vein, at 
fixed planet mass, the spread in the migration timescale 
should increase with decreasing /3. This trend does not 
clearly show up in Fig. [4j most probably because our 
minimum and maximum values of /3 only differ by a fac- 
tor of 2. Some further insight into the effect of stochastic 
kicks is given in Section [3.2.21 We finally point out that 
the disc's density and temperature profiles do not signifi- 
cantly decrease over the planets migration timescale (for 
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Figure 4. Time evolution of the orbital separation of a Saturn-mass planet (left column), a Jupiter-mass planet (middle column), 
and of a 5 Jupiter-mass planet (right column) in a gravitoturbulent disc with /3 = 15 (top row), = 20 (middle row), and f} = 30 
(bottom row). Time is displayed in years (bottom x-axis), and in orbital periods at 100 AU (top x-axis). Results are obtained for 
eight evenly-spaced planet's azimuths at restart (azimuth increases with wavelength's color). 



the longest simulation with the Saturn-mass planet, the 
disc's mass decreases by at most 10%). 

3.2.2 Stochastic kicks 

The results of S ect ion \3 . 2 . 1 1 show that planets formed by 
gravitational instability should rapidly migrate inwards. 
Even for the Saturn-mass planet, whose mass is smaller 
than the estimated initial clump's mass (see first para- 
graph of Section I3.2|l , inward migration occurs in typi- 
cally less than 10 orbits. As illustrated in Fig. 2] planet 
migration in gravitoturbulent discs is not a smooth pro- 



cess. Planets may experience large kicks, either inwards or 
outwards, depending on the local density perturbations 
they encounter. An example of outward kick is illustrated 
in Fig. [5] which displays the gas density obtained for 
the restart simulation with /3 = 30 and the Saturn-mass 
planet that experiences a significant episode of outward 
migration in the bottom-left panel of Fig. [4] (red curve) . 
Time increases by 0.1 orbital periods at 100 AU from left 
to right. In each panel, the planet's location is indicated 
with a filled black circle. 

We first comment that, because the planet's mass 
is low enough (q/hp ~ 0.3), the wakes it generates are 
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Figure 5. Contours of the gas surface density obtained for the restart simulation with /3 = 30, and the Saturn- mass planet with 
restart azimuth equal to 7tt/4. This simulation corresponds to the red curve in the bottom- left panel of Fig. [4] where a significant 
episode of outward migration occurs from 3 to 6 orbits after restart. Orbital separation is shown on the x-axis, and azimuth on 
the y-axis. In all panels, the planet's location is depicted with a filled black circle. 



invisible, their density contrast being much weaker than 
that of the turbulent perturbations. In the case depicted 
in Fig. [5] the planet finds itself with an underdense re- 
gion behind it, and a (comparatively) overdense region in 
front of it. This density contrast yields a large, positive 
coorbital corotation torque that is responsible for the vig- 
orous outward kick experienced by the planet. This kick 
increases the planet's orbital separation from ~ 75 AU to 
~ 95 AU in only 0.2 orbital periods at 100 AU (that is, 
in 200 years). Similarly, inward kicks are obtained as the 
planet experiences a negative coorbital corotation torque 
induced by an underdense region ahead of it, and a (com- 
paratively) overdense region behind it. As is shown for 
instance in the bottom-left panel of Fig. [4] kicks excite 
planets epicyclic motion, which is rapidly damped as the 
underlying disc-planet tidal torque takes over. 

The stochastic kicks induced by gravitoturbulence 
are reminiscent of type III runaway migration. This mi- 
gration regime relies on the formation of a coorbital mass 
deficit in the planet's horseshoe region, which takes the 
form of an asymme tric density structure behin d and in 
front of the planet (|Masset fc Papaloizoul 120031 ). In mas- 
sive (Q > 1) laminar discs with moderate viscosity (a 
is a few xlO -3 ), runaway migration is particularly rele - 
vant for planets with q ~ hp l|Masset fc Papaloizoull2003l ) , 
like the planets considered in our study. However, in our 
model, the disc's turbulent viscosity (atot is a few per- 
cent) is too large to allow planets to open a partial dip 
around their orbit before they reach the disc's inner edge. 
Thus, planets do not build up a coorbital mass deficit. 
But, through the density fluctuations it triggers, gravito- 
turbulence may provide an effective mass deficit on each 
side of the planet, and induce an effective type III migra- 
tion. However, this effective mass deficit is not coorbital 
with the planet (as can be seen in Fig. [5} , and therefore 
cannot lead to a runaway. Consequently, the relative fast 
inward migration that we obtain cannot be attributed to 



type III runaway migration. Also, although the planets 
we consider are massive, their short formation and migra- 
tion timescales compared to their gap-opening timescale 
suggests they are subject to type I-like migration instead 
of type II, perturbed by stochastic kicks. This point will 
be clarified in Section \3. 2. 31 

Stochastic kicks provide a source of very fast relative 
motion between the planet and the disc. Fig. [6] displays 
the migration speed in units of the local sound speed 
for the three runs with ft = 30, and zero azimuth at 
restart. The migration rate can be a significant fraction 
of the local sound speed, and short episodes of supersonic 
motion are even obtained. We note that the maximum 
value for the migration speed to sound speed ratio tends 
to increase with decreasing planet mass, as less massive 
planets are more sensitive to stochastic kicks. 

3.2.3 Comparison to laminar disc models 

Our results of simulations indicate that the migration 
of planets embedded in gravitoturbulent discs is driven 
by both a vigorous (negative) tidal torque, and a (posi- 
tive or negative) stochastic torque. As the planet mass 
increases, the former prevails over the latter. In this 
section, we compare the results of our turbulent disc 
model with f3 = 30 to those of an equivalent laminar 
disc model, including the viscous force in the momentum 
equation. A constant alpha viscosity is used, approxi- 
mately equal to 1.3%. The comparison is done for the 
same three planet masses as previously considered. As 
the initial conditions for the laminar calculations, we take 
surface density and temperature profiles that fit those 
of the gravitoturbulent run with /3 — 30 at 30 orbits: 
£ « 15 g cnr 2 x (r/100 AU)~ 3/2 , and T « 25 K. 

Two series of laminar calculations were carried out: 
(i) one with an isothermal equation of state, where the 
imposed (flat) temperature profile remains constant in 



© 2011 RAS, MNRAS 000.fTHT3l 



8 



C. Baruteau et al. 



Time (orbits after restart) 
2 3 4 




1000 



2000 3000 
Time [years] 



4000 5000 



100 



3 
< 



o 



CO 
Q. 


in 



Time (orbits after restart) 
10 15 




5.0-10 



1.0-10 4 1.5-10 4 
Time [years] 



2.0-1 4 2.5-1 4 



Figure 6. Time variation of the (absolute value of the) migra- 
tion speed, in units of the local sound speed, obtained for the 
restart simulations with /3 = 30 and zero azimuth at restart. 
They correspond to the results shown with black curves in the 
bottom row of Fig. [4] The migration speed is measured as the 
time derivative of the planet's orbital separation. 

time, and (ii) one with an energy equation with both 
viscous heating (a ~ 1.3%) and /3— cooling (/3 = 30). For 
the latter series, the aforementionec0 value for a leads to 
an initial thermal balance between viscous heating and 
/?— cooling, see Eq. @. 

All laminar calculations include the axisymmetric 
component of the disc's self-gravity. This is meant to 
avoid a large mismatch between the disc's and the 
planet's angular velocities due to the disc gravity, which 
would yield a spurious shift inwards of all Lindblad reso- 
nances, and ther efore an artificial accelerat ion of the in- 
ward migration ()Baruteau fc Massetll2008bl ). Also, in all 
laminar runs, the calculation of the force exerted by the 
disc on the planet excl udes the planet's circumplanetary 
disc (|Crida et al.ll2009h . All other parameters are other- 
wise the same as in the gravitoturbulent calculations. 

The results with the laminar disc model and an en- 
ergy equation are depicted in Fig. [7] as solid curves. It 
takes about 5 orbits for the 5 Jupiter-mass planet to mi- 
grate from 100 AU to 50 AU, and about 7 orbits for the 
Jupiter- mass planet. This corresponds to migration rates 
~ half those of the gravitoturbulent calculations (see bot- 
tom row of Fig. [4|. Interestingly, the difference in mi- 
gration rates is larger for the Saturn-mass planet, which 
drifts from 100 to 60 AU in 20 — 25 orbits in the laminar 
run, and (on average) in only 5 — 6 orbits in the gravito- 
turbulent run. We note that the migration timescales of 
the Saturn- and Jupiter-mass planets differ by a factor 
approximately equal to their mass ratio, as expected in 
the type I migration regime, wherein the tidal torque en- 
compasses the differential Lindblad torque and the horse- 
shoe drag (|Paardekooper Sz Papaloizoul 120091 '). We com- 
ment that it is the large disc aspect ratio (h ~ 0.1) 

3 As in the case with = 20, this value of a, given by viscous 
transport theory, is slightly smaller than what we measured 
in the gravitoturbulent disc model ((«tot) ~ 1.8%). We have 
checked that a higher value of a does not change our results 
of laminar calculations, presented below. 



Figure 7. Time evolution of the planets' orbital separation 
obtained for the laminar disc model with initial profiles of 
surface density and temperature similar to those of the gravi- 
toturbulent run with j3 = 30 before restart. Solid curves show 
the results with inclusion of an energy equation (with vis- 
cous heating and /3-cooling), and dashed curves those with an 
isothermal equation of state. 

that makes the Saturn-mass planet (q/h 3 ~ 0.3) and the 
Jupiter-mass planet (q/h 3 ~ 1) relevant to type I migra- 
tion. The 5 Jupiter-mass planet (q/h 3 ~ 5 ^> 1) is, how- 
ever, more prone to non-linear effects (e.g.. iMasset et ail 
2006). This is presumably why the migration timescales 
of the Jupiter- and the 5 Jupiter-mass planets do not dif- 
fer by a factor of their mass ratio. It may seem surprising 
that the 5 Jupiter-mass planet does not open a gap and 
migrates in a very short timescale. The gap-opening cri- 
terion for a planet on a fixe d circular orbit in a laminar 
viscous disc takes the form (|Crida et al.ll200(3 ) 

where a and h are to be evaluated at the planet's or- 
bital separation. Although it is unclear to what extent 
the gap-opening criterion in Eq. © may be applied to 
planets migrating in gravitoturbulent discs, it provides 
some insight into the fact the planets in our disc models 
are not expected to clear a gap. Notwithstanding that 
the large values of h and a imply the 5 Jupiter-mass 
planet does not quite satisfy the above gap-opening cri- 
terion (the left-hand side of Eq. (J3]) typically ranges from 
2 to 4 in the disc inner parts), we stress that the vig- 
orous tidal torque exerted by the disc, which is partly 
due to the large masses of the disc and the planet, con- 
spires to make the planet migrate in a timescale shorter 
than the timescale required to open a gap (a few libra- 
tion timescales, which would correspond here to a few 
tens of orbits). In other words, the 5 Jupiter-mass planet 
migrates very fast because it has no time to open a gap. 

The results of laminar simulations with an isother- 
mal equation of state are overplotted as dashed curves 
in Fig. [7] The migration rates obtained in this se- 
ries are smaller than with inclusion of an energy 
equation. The difference arises from a large, negative 
horseshoe drag obtained with an energy equation, and 
driven by a positive entropy gradient in our disc model 
l|Baruteau fc Massetll2008al ; iPaardekooper fc Papaloizoul 
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Figure 8. Time variation of the mass enclosed in the planet's 
Hill radius in units of the planet mass, obtained for the restart 
simulation with = 30 and and zero azimuth at restart. The 
mass ratio is displayed before the planet crosses the inner edge 
of the computational domain. 

l2008l : lMasset fc Casolill2009l ; iPaardekooper et al ] |20ld . in 

our case, the radial profile of the disc entropy scales in 
proportion to the orbital separation). 

Further insight into the impact of the horseshoe drag 
may be obtained by comparing the migration timescales 
of our laminar calculations to those predicted by type 
I migration theory. Note that the predicted migration 
timescales are inferred from the torque exerted on a 
planet held on a fixed circular orbit. In the isothermal 
series, where E tx r _3//2 and T is uniform, the total 
torque reduces to the differential L i ndbla d torque. Us- 
ing Eq. (14) of Paardcko oper et al.l (|2010h . we find mi- 
gration timescales of about 75 and 20 orbits for the 
Saturn- and Jupiter-mass planets, respectively. This is 
in good agreement (to within ~ 25%) with our find- 
ings (dashed curves in Fig. [7}. In the series with the 
energy equation, the total torque comprises the differ- 
ential Lindblad torque and the horseshoe drag. Despite 
the large value of a in our disc model, the horseshoe 
drag is found to be fully unsaturated, as the diffusion 
timescale across the planet's horse shoe region remains 
larger than the U-turn timescale (IBaruteau fc Massetl 
l2008al ; [Masset fc Casohll2010l: IPaardekooper et alJj201lT) . 
Using Eqs. (14) and (45) of IPaardekooper et alj (|201Ch . 
we now find migration timescales of about 35 and 10 
orbits for the Saturn- and Jupiter-mass planets, respec- 
tively. This is shorter than what we find (see solid curves 
in Fig. [7), although still in decent agreement. The differ- 
ence presumably arises from the large velocity difference 
between the disc and the planet, which sets asymmetric 
horseshoe stre amlines around the planet's orbital radius 
([Massed I2OO2T ), an effect that is not taken into account 
i n the horseshoe drag expression of iPaardekooper et al.l 
(|2010h . 

We further comment on having a positive entropy 
gradient in our gravitoturbulent disc model. Let us as- 
sume the disc profiles of surface density and temperature 
scale as r~ a and , respectively. The disc entropy s, 
which we define as s = pE -7 , where the disc pressure 
p verifies the ideal gas equation of state, scales as r -5 , 



with £ = w — (7 — l)cr. As the disc evolves towards a 
steady state with approximately uniform Toomre-Q pro- 
file, £ ~ —3 + a(3 — 7). Steady-state gravitoturbulent 
discs may thus have a globally decreasing profile of en- 
tropy only for surface density profiles steeper than r~ 9 ^ 4 , 
assuming 7 = 5/3. Furthermore, for 7 = 5/3 and our 
fiducial softening length 0.3H(r p ), the fully unsatu rated 
horseshoe drag Fh s reads ( Paardcko oper et al.ll2010l . their 
Eq. 45) 

7lWr ~ -20.5 + 8.6a, (4) 

where To = (q/h p ) 2 T, p rpQ,p, and all quantities with the 
subscript p are to be evaluated at the planet's location. 
Eq. Q shows that, in discs with a uniform Toomre-Q 
parameter, the horseshoe drag is negative whenever the 
surface density profile is shallower than r~ 2A , assuming 
7 = 5/3. The dif ferential Lindblad torque T l can be sim- 
ilarly recast as l|Paardekooper et al. I I2010I . their Eq. 14) 

7 r L /r ~ 3.2 - 4.0a, (5) 
and the total torque r to t as 

7r t ot/r = -17.3 + 4.6a. (6) 

Using Eq. JBJl, the type I migration timescale, 77 = 
rpfipMp/2r t ot, reads 

£«« »<"-.>-'■*. (4)" (&)"". m 

where T or b stands for the orbital period at the planet's 
orbital radius. The estimate in Eq. (0 applies to planets 
with q ~ hp, and is approximate, since it is based on the 
torque expression for a planet on a fixed orbit. As we have 
already stressed above, the fast relative motion between 
the disc and planet is likely to affect this estimate. Nev- 
ertheless, it indicates that the type I migration of planets 
formed by fragmentation (q ~ hp, Q p ~ 1.5) in discs with 
moderate aspect ratio (hp ~ 0.1) should typically occur 
in a few tens of orbits at 100 AU, that is a few 10 4 years. 

Still, the migration timescales in our gravitoturbu- 
lent simulations are a factor of a few shorter than ob- 
tained in laminar disc models. In addition to the stochas- 
tic kicks due to turbulence, the following two mechanisms 
may account for the larger migration rates in the gravi- 
toturbulent disc model: 

- The non-axisymmetric component of the disc's self- 
gravity shifts L indblad resonances tow ards the planet's 
orbital radius l|Pierens fc Hurel 120051 ). For discs close 
to marginal stability, and for our softening parameter 
(e/H ~ 0.3), this shift increases t he differential Lind- 
blad to rque by typically a factor of 2 (|Baruteau fc Ma sset 
2008b, their Figure 10). Note however that this mech- 
anism has been investigated in laminar disc models. It 
is therefore not straightforward whether it should be- 
have identically with turbulence. Still, numerical simu- 
lations indicate that, on average, the Lindblad torque 
behaves similarly in la minar and turbulent disc models 
(jBaruteau fc Linll2010l . where wave-like turbulence is gen- 
erated by stochastic forcing of the disc). 
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Figure 9. Migration of a Jupiter-mass planet in a disc with cooling parameter f) = 30. In the left panel, the planet mass is 
gradually ramped up over one orbit, while the planet is held on a fixed circular orbit. After this first stage (not shown in this 
panel), the planet is released in the disc. The middle panel shows results at higher resolution (1024 X 2048), and the right panel 



shows results with increased self-gravitating softening length (e s 



3 X 10 r). All three panels should be compared to the 



mid-bottom panel in Fig- [4] which displays the results of our fiducial setup (no gradual increase of the planet mass, grid resolution 
is 512 X 1536, and s sg = 3 X 10" 4 r). 



- In fully self-gravitating discs, the effective mass that 
dictates the strength of the tidal torque and thus the mi- 
gration rate not only includes the planet mass, but also 
the mass of the gas surrounding the planet's location. For 
massive planets (q > hp), the latter corresponds to the 
planet's circumplan etary disc, which is a fraction of the 
planet's Hill sphere (|Crida et al.ll2009l ). It is unclear what 
it corresponds to for smaller planets; we speculate that it 
might be the planet's Bondi sphere. For illustration pur- 
poses, the time evolution of the mass inside the planet's 
Hill radius is shown in Fig. [8] for three gravitoturbulent 
runs with f3 = 30. It is comparable to, albeit somewhat 
larger than the planet mass. It progressively grows as 
the planet moves inwards, despite the decrease in the 
planet's Hill radius with decreasing orbital separation. It 
is possible that the short cooling timescales used in our 
simple cooling prescription (a few dynamical timescales, 
depending on 0) implies that the gas in the planet's Hill 
radius can contract more rapidly than it would do with 
a more realistic cooling function. Again, notwithstanding 
the uncertain dynamical role of the planet's Hill sphere in 
turbulent disc models, it is likely that the increased effec- 
tive mass is responsible for an increase in the migration 
rate by a factor of a few compared to laminar calculations 
(which only include the axisymmetric component of the 
disc's self-gravity). We comment that both the increase 
in the effective planet mass, and the decrease in the disc 
aspect ratio with radius (in discs with fiat temperature 
profile, h oc r 1 ^ 2 ), make planets less and less sensitive to 
stochastic kicks as they migrate inwards. This explains 
why in Fig. [S] the amplitude of the stochastic kicks de- 
creases with decreasing orbital separation. This also im- 
plies that, if a planet undergoes a vigorous inward kick, it 
will be increasingly harder to reverse the trend. However, 
should a planet experience a large outward kick, it would 
take longer for the tidal torque to take over and lead to 
inward migration again. 



3.2.4 Impact of numerics 

We briefly assess in this section the impact of numer- 
ics on our results of gravitoturbulent calculations. We 
have performed an additional suite of simulations with 
the Jupiter-mass planet and ft = 30, with the planet's 
mass first smoothly ramped up over one orbit, while the 
planet is held on a Keplerian orbit. After this, the planet 
is released in the disc. This procedure is meant to mimic 
the clump formation, which typically occurs over a dy- 
namical timescale, much shorter than the gap-opening 
timescale (see Section 13 . 2 - 3 p . The time evolution of the 
planet's orbital separation is depicted in the left panel of 
Fig. [9] It shows the very same trend as obtained with- 
out smoothly increasing the planet's mass (see the mid- 
bottom panel of Fig. [3} . 

To assess the influence of resolution on our simula- 
tion results, we have varied the grid's resolution, and the 
softening length £ sg over which the self-gravitating ac- 
celeration is smoothed. We have run two additional disc 
models with /3 = 30, one with an increased grid resolu- 
tion (1024 x 2048), and one with a larger e sg (3 x 10~ 2 r). 
In the latter 0.3 times the initial pressure 

scale height, and it can be regarded as a crude treat- 
ment for the effects of the disc's finite thickness. The 
value of e S g at r = r p equals the softening length of the 
planet's gravitational potential. We have checked that be- 
fore the introduction of the planet at 30 orbits the disc 
properties (density and temperature profiles) in both ad- 
ditional models are in close agreement with those of the 
corresponding run at our fiducial resolution and soften- 
ing length value. We restarted these simulations with a 
Jupiter-mass planet. Again, in each case, a series of eight 
simulations was performed with evenly-spaced values for 
the planet's azimuth at restart. Results with increased 
grid resolution, and increased e sg are depicted in the 
middle and right panels of Fig. [9] respectively. Overall, 
the amplitude of the stochastic fluctuations is somewhat 
larger in the high-resolution simulations. Still, except in 
one run that features a rather long episode of outward 
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Figure 10. Time evolution of the planets' orbital separation for the disc model with ft = 15 restarted at 85 orbits, when the disc 
aspect ratio is ss 6%. From left to right, results are shown for a 20 Earth-mass planet (q = 6 X 10 — 5 , q/h 3 ~ 0.3), a ~ Saturn-mass 
planet (q = 2 X 10 _4 ,g//i 3 ~ 1-5), and a Jupiter-mass planet (q = 10 -3 , q/h 3 ~ 5). As in the previous figures, results are obtained 
for eight evenly-spaced planet's azimuths at restart (increasing with increasing wavelength's color). 



migration, the migration timescale of the high-resolution 
simulations ranges from about 3 to 9 orbits, which is 
broadly consistent with the results at our fiducial res- 
olution. Furthermore, increasing e sg from 3 x 10 -4 r to 
3 x 1CP 2 r does not significantly change the amplitude of 
the stochastic kicks experienced by the planet, and the 
migration timescales obtained with both values of e sg are 
in reasonable agreement. The results of both additional 
disc models indicate that the planet's orbital evolution 
is essentially unaltered by density structures arising from 
gravitoturbulence with typical size < 0.3H. 



3.3 Results at lower disc temperature 

As shown in Section T3. 21 the migration of planets formed 
by gravitational instability results from a competition 
between the background disc-planet tidal torque, and 
the stochastic torque due to turbulence. At least in the 
early stages of their evolution, such planets should have 
a planet-to-primary mass ratio q ~ h p , and be subject 
to type I- like migration (see Section f3. 2. 3jl . As far as the 
tidal torque is concerned, Eq. shows that the migra- 
tion timescale is essentially controlled by three dimen- 
sionless parameters: q/hp, which should be ~ 1 in the 
early stages of the planet evolution, Q p , which is > 1.5 
in a steady state, and h p . The disc model that we have 
considered so far has h p ~ 0.1, which leads to rapidly 
formed Jupiter-mass planets. However, it is possible that 
the disc initially has a lower aspect ratio, thereby forming 
less massive planets. The latter may still become Jupiter- 
sized objects if mass growth is efficient. 

To evaluate the impact of a lower disc temperature 
on our results, we restarted the simulation of Section [3T] 
with /3 = 15 at 85 orbits. At this time, the disc aspect 
ratio at the planet's restart location (100 AU) has de- 
creased from 0.1 to about 0.06, and the disc's mass from 
0.2M* to 0.15M*. We choose again three planet masses 
corresponding to a planet-to-primary ratio q around h p : 
a 20 Earth-mass planet (q = 6 X 10' b ,q/h 3 p ~ 0.3), a 
~ Saturn-mass planet (q = 2 x 10~ 4 , q/h p ~ 1), and a 
Jupiter-mass planet (q = 10 -3 , q/h p ~ 5). As previously, 



we take eight evenly-spaced planet's azimuths at restart. 
The time evolution of the planets' orbital separation is 
shown in Fig. 1101 

Compared to the results with h p =0.1, planets with 
the same q/h p are now more sensitive to stochastic fluc- 
tuations, the tidal torque becoming less vigorous. This 
can be seen from Eq. (0, which indicates that at fixed 
values of q/h p and Q p , the migration timescale due to 
the tidal torque scales with h p 2 . The mean migration 
timescales that we obtain for h p ~ 0.06 (about 6, 12, and 
50 orbits with decreasing q/h p ) are typically a factor of 
3-5 longer than in the corresponding runs with h p ~ 0.1, 
which is roughly consistent with the aforementioned scal- 
ing Tjr oc h p 2 . 



4 CONCLUDING REMARKS 

Planets observed at large orbital separation from their 
host star (typically > 50 AU) are thought to be potential 
candidates for the formation scenario based on the gravi- 
tational instability of massive discs. A number of studies 
have shown that the formation of clumps by fragmenta- 



tion is possible at such separations (e.g., Rankovll200 



iMatzner fc Levin|[2005l ; IStamatellos fe Whitworthll200; 
The orbital evolution of clumps formed by gravitational 
instability has be e n investiga t ed in a few studies (e.g. . 



Maver et alj 12002 



Bossll2005l; IVorobvov fc Basul l2010al : 



Cha fc Navakshinl l201oi iMachida et all l201ll ). Clumps 



are often observed to migrate inwards on short timescales, 
although in some disc models the formation of several 
clumps and/or the use of different cooling functions has 
a significant impact on migration, which sometimes ap- 
pears to be suppressed l|Bossll2005l ). 

In this paper, we have studied the interaction be- 
tween a single planet and the gravitoturbulent disc it is 
embedded in, following the assumption that the planet 
has formed by gravitational instability. We have per- 
formed 2D hydrodynamical simulations with a simplified 
prescription for the disc's cooling. Three disc models with 
different cooling timescales have been considered, giving 
rise to different levels of gravitoturbulence. After having 
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set up a quasi-steady state gravitoturbulent disc, we have 
restarted our simulations including a single planet with a 
range of masses approximately equal to the expected ini- 
tial mass of clumps formed by fragmentation (> h 3 M t , 
wit h h the disc aspec t ratio at the fragmenting location; 
see lBolev etalll2010h . 

We find that a planet interacting with its nascent 
gravitoturbulent disc migrates inwards on very short 
timescales, despite the stochastic kicks due to turbulent 
density fluctuations. In a disc with aspect ratio ~ 0.1, 
it takes less than 10 orbits at 100 AU (that is, about 
10 4 years) for a Jupiter-mass planet to migrate from 
100 AU to 20 AU (the inner edge of our computational 
grid). Planets less massive than the expected fragment- 
ing mass are more sensitive to stochastic kicks, but also 
migrate inwards relatively fast. Their formation and mi- 
gration timescales being shorter than their gap-opening 
timescale (partly because of the disc's vigorous turbu- 
lence), planets formed by fragmentation do not open a dip 
or a gap around their orbit, and are therefore not subject 
to the type II and runaway type III migration regimes. 
Turbulent density fluctuations may provide an effective 
mass deficit in the planet's horseshoe region that kicks 
the planet either inwards or outwards (see Section [3]2]2}. 
These kicks can be seen as an effective, temporary type 
III migration coming on top of the inward migration due 
to the background disc-planet tidal torque. 

The comparison of our results with those of laminar 
disc models shows that the averaged tidal torque driv- 
ing the net inward migration in gravitoturbulent discs 
is essentially the one that would lead to type I migra- 
tion in the absence of turbulence (see Section T3. 2. 3 [I . Part 
of the reason as to why planets rapidly migrate inwards 
in our model is due to the set up of a positive entropy 
gradient in a steady state, which yields a large nega- 
tive horseshoe drag, adding up to the negative differ- 
ential Lindblad torque. We argue that in gravitoturbu- 
lent discs with uniform Toomre-Q profiles and cooling- 
to-orbital timescale ratios (that is, with a uniform alpha 
parameter), the horseshoe drag should be negative, un- 
less the surface density decreases very steeply (see Eq. [4j . 
Although this possibility cannot be excluded, we believe 
that fast inward migration should be a generic expecta- 
tion for planets formed by gravitational instability. 

The main conclusion of this present study is that 
a single planet formed by gravitational instability should 
migrate inwards on a timescale much shorter than the ex- 
pected lifetime of protoplanetary discs. However, it is not 
possible at this stage to say that massive planets at large 
separation from their host star, such as the HR 8799 sys- 
tem, could not have formed by gravitational instability. 
The simulations presented here contain a number of sim- 
plifications that need to be addressed. We have consid- 
ered a simplified cooling prescription, causing the entire 
disc to be gravitoturbulent. In reality, we would expect 
the inner parts of a disc to be too hot to be gravitation- 
ally unstable, and other sources of turbulence such as 
the magnetorotational instability could prevail (see e.g., 
IZhu et al.ll2009l ; lRice fc Armitagell2009l ; IClarkell2009l , who 
studied the radial dependence of the alpha viscosity pa- 
rameter in self-gravitating discs with realistic cooling). It 
is thus possible that the rapid type I migration of planets 



formed by gravitational instability gets slowed down in 
the disc inner parts, resulting in gap formation. A more 
realistic cooling prescription will also impact the rate at 
which planets can contract, with various outcomes for 
the evoluti on of the planet mass over its orbital migra- 
tion (e.g.. lNavakshiij|2010h . Also, we have focused on 
the orbital evolution of a single planet assumed to have 
formed by fragmentation, whereas the formation of multi- 
ple clumps is a likely outcome of gravitational instability. 
It is fully expected that the interaction between clumps 
may significantly impact their migration. We will address 
these issues in future work. 



APPENDIX A: CALCULATION OF THE 
ALPHA PARAMETERS IN A 
TWO-DIMENSIONAL DISC 

When mass is confined to a razor-thin disc, the cal- 
culation of the alpha parameters associated with the 
Reynolds and gravitational stresses may be inferred from 
the general three-dimensional case by substituting the 
volume mass density p with E<5(z), where E is the sur- 
face mass density of the razor-thin disc, S is Dirac's delta 
functi on, and z denotes the vertica l coordinate. O ne ob- 
tains l|Lvnden-Bell fc Kalnaislll972l ; lGammidl200lh 



CtRcy 



2 (T,5v r 5v v ) 

3 (Ecf) 1 



(Al) 



where (.) stands for the azimuthal average, 5v r = v r — 
(v T ), 5v v = v v — (Vtp), and 

2 {JZ (^G)- 1 g r g p dz) 
Qss= 3 (Ed) ' (A2) 

In both previous equations, v r , v v and c s are to be eval- 
uated in the razor-thin disc (z = 0) and are therefore 
directly computed in a 2D simulation. However, the ver- 
tical dependence of g r and g v , the radial and azimuthal 
components of the self-gravitating acceleration, needs to 
be determined as the gravitational field outside the disc 
contributes to the shear stress despite mass being con- 
fined to the disc. The potential $ from which g r and g v 
can be derived is given by 



4> 



vA*2 + r r - 



■ 2rr' cos(<^ — ip') + z 2 + e| g 



(A3) 



where T> is the domain where E does not vanish, and 
£ sg is a small softening length required to avoid numeri- 
cal divergences when z = 0. From the surface density of 
our 2D simulations, g r (r, ip, z) and g v (r, ip, z) can be con- 
veniently calculated by fast Fourier transforms provided 
that z is taken proportional to r, condition required for 
g r and g v to read as convolution products (another nec- 
essary condition is that a grid with logarithmic radial 
spacing is used). The same condition applies to £ sg for 
the calculation of g r and g v at z = 0. Writing e sg = Br 
and z — nr, with B and r\ constants, the expressions for 
q r (r, lp, z) and g^(r, <P, z ) are gi ven by Eqs. (Al) and (A3) 
of iBaruteau fc Massed (20083) with B' 2 -> B 2 + r/ 2 . The 
vertical integral in Eq. (|A2I) is then numerically calcu- 
lated from the values of g r and g v obtained with a range 
of values for r\. 
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